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ABSTRACT 

We use time-varying models of the coupled evolution of the HI, H2 gas phases and stars 
in galaxy-sized numerical simulations to: a) test for the emergence of the Kennicutt-Schmidt 
(K-S) and the H2-pressure relation, b) explore a realistic H2-regulated star formation recipe 
which brings forth a neglected and potentially significant SF-regulating factor, and c) go be- 
yond typical galactic environments (for which these galactic empirical relations are deduced) to 
explore the early evolution of very gas-rich galaxies. In this work we model low mass galaxies 
C^baryon < 10 9 M Q ), while incorporating an independent treatment of CO formation and de- 
struction, the most important tracer molecule of H2 in galaxies, along with that for the H2 gas 
itself. We find that both the K-S and the ^-pressure empirical relations can robustly emerge 
in galaxies after a dynamic equilibrium sets in between the various ISM states, the stellar com- 
ponent and its feedback (T > 1 Gyr). The only significant dependence of these relations seems 
to be for the CO-derived (and thus directly observable) ones, which show a strong dependance 
on the ISM metallicity. The H2-regulated star formation recipe successfully reproduces the 
morphological and quantitative aspects of previous numerical models while doing away with 
the star formation efficiency parameter. Most of the HI — > H2 mass exchange is found taking 
place under highly non-equilibrium conditions necessitating a time-dependent treatment even 
in typical ISM environments. Our dynamic models indicate that the CO molecule can be a poor, 
non-linear, H2 gas tracer. 

Finally, for early evolutionary stages (T < 0.4 Gyr) we find significant and systematic devi- 
ations of the true star formation from that expected from the K-S relation, which are especially 
pronounced and prolonged for metal-poor systems. The largest such deviations occur for the 
very gas-rich galaxies, where deviations of a factor ~ 3 - 4 in global star formation rate can 
take place with respect to those expected from the CO-derived K-S relation. This is particu- 
larly important since gas rich systems at high redshifts could appear as having unusually high 
star-formation rates with respect to their CO-bright H2 gas reservoirs. This points to a possi- 
bly serious deficiency of K-S relations as elements of the sub-grid physics of star formation in 
simulations of structure formation in the Early Universe. 

Subject headings: galaxies: numerical simulations - galaxies: spirals - galaxies: star formation 
- ISM: molecular gas - ISM: atomic gas - molecules: H2,CO 
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1. Introduction 



In spite of the fact that the general character of the cycle through which galaxies convert their ISM 
to stars has been known for a long time, it has proven to be remarkably difficult to formulate a predictive 
theoretical framework for this process. Indee d, we know that throughout most of the Univers e star formation 
takes place in molecular gas complexes (e.g. lSolomon and Vanden Boual2005uOmondl2007r) . The HI — > H2 
phase transition is conditioned by a combination of sufficiently high HI column densities and pressures, 
and consequently star formation tends to concentrate in high density regions, e.g. in the central parts of 
galaxies, in spiral arms or high-pressure concentrations of gas formed by bulk gas motions or swept up 
by the shocks from supernovae and stellar winds of OB associations. The link between molecular gas and 
star formation is so tight that it has even been used to infer the distribution of the former by that of th e 
latter, when the CO-H2 conversion factor was still considered v ery uncertain ( Rana and WilMnsorJ ll986). 
In the Galaxy, where this link is best studied (e.g. lBlitall997l and references therein), the latest results 
confirm st ar formation always ta king place in CO-bright molecular clouds, even at very large galactocentric 
distances (IKobayashi et al. In other ga laxies this tight associ ation has been verified in all cases where 

sufficient angular resolution is available (e.g. lWong and Blital2002r) . Thus it is fair to say that ^formation 
is a necessary prerequisite for star formation in galaxies, and incorporating it in galaxy-sized numerical 
simulations of gas and stars is the single most important step currently missing from a realistic rendering of 
star formation in such models. 

Following the early and widespread observational evidence establishing the H2-(star formation) link. 



the inclusion of the H2 gas phase in numerical models of galaxies has occured only recently (IPelupessy et al. 



2006: iDobbs et all 12006 : iRobertson and Kravtsovl 120081) . and the refinement of these models is an area of 
ongoing research. This is mainly due to the difficulty of tracking the dynamic and thermodyna mic evolution 
of H? and its precursor phase, the Cold Neutral Medium HI (n ~ 5 - 100 cm -3 , T^ ~ 60-200K. IWolfire et al. 



l2003n in galaxy-sized numerical models and due to the strong H2 self-shielding complicating radiative trans- 
fer models of its far-UV radiation-induced destruction. The first problem has prevented most efforts from 
properly tracking the HI — > H? phase transition in galaxies wit hout resorting to simplifying steady-state 
solutions ( Hidaka and Sofuell2002 : Robertson and Kravtsov 2008) (suitable only for quiescent galactic envi- 
ronments), or to semi-empirical multiphase models (e.g. lSemelin and Combesll2002l) with limited predictive 
value. The second problem confounds even numerical simulations tracking the HI — ► H2 phase transition in 
individual gas clouds where local approximations of the self-shielding HI/H2 volume (necessary for numer- 
ically managea ble solutions) can make the H2 gas mass fraction a strong function of the chosen numerical 
resolution (e.g. iGlover and Mac Lowll2007l) . Finally a secondary, yet important problem of such single gas 
cloud simulations is posed by the constant boundary conditions assumed during their evolution, which are 
an unlikely setting for real gas clouds immersed in the ISM environment of a galaxy. In such environments 
cloud boundary conditions that powerfully influence the HI — > H2 phase transition, such as the ambient FUV 
radiation field and pressure, change on timescales comparable or shorter than "intern al" cloud dynamic and 
thermodynamic timescales, especially in vigorously star forming environments (e.g. 
Wolfire et al.ll2003l : IPelupessv etaDbood) . 



Parravano et al 



2003; 



Despite the aforementioned difficulties the incorporation of the H2 <-> HI gas phase interplay, and its 



-3- 



strong role as star formation regulator, in numerical models holds the promise of large improvements in their 
handling of galaxy evolution, and the possibility of unveiling new, hitherto neglected, aspects of star for- 
mation feedback o n the ISM. In this pap er we apply our numerical models for the coupled evolution of gas 
(HI, H2) and stars dPelupessy et al.l l2006) to the evolution of low mass galaxies (Mbaryon < 1O 9 M ) in order 
to explore two new key directions, the first of which is the emergence of two important empirical relations 
found for galaxies in the local Universe: the Kennicutt-Schmidt (K-S) and the ^-pressure relation. Sec- 
ondly we will make a investigation of very gas-rich systems (more typical of the Early Universe), and check 
whether the aforementioned relations remain valid during their evolution. The latter is of crucial importance 
given the prominant role such empirical relations are given in describing the sub-grid star-formation/gas 
interplay in cosmological simulations of galaxy evolution (where the resolution limitations imposed by the 
simulation of large volumes preclude a detailed description of star formation). Finally along with our orig- 
inal time-dependent treatment of the HI — > H2 phase transition we also include the CO molecule, allowing 
direct comparisons to the observed H2 gas distributions, and a new independent investigation of the CO-H2 
relation within the dynamical setting of an evolving galaxy. 

The structure of our work is as follows: in section [2] we present relevant features of the model, show 
semi-analytical predictions for the ^-pressure relation, and formulate the extention of the model so that 
includes CO, in section [3] we present our detailed numerical simulations, and investigate the K-S, H2- 
pressure and CO-H2 empirical relations. In section|4]we investigate and discuss the validity of the important 
K-S relation during early galaxy evolution stages, and for very gas-rich systems. We then summarize our 
conclusions in section |5] 



2. Incorporating the molecular gas phase: a dynamical approach 



The HI — ► H2 phase transition in galaxies, as catalyzed by dust grains, have been studied extensively 
ever sinc e the strong self-shielding nature of H2 in its dissociation by far-U V (FUV) photons ha s been rec 



ognized ( Spitzer and Jenkins 



1975 



Sava ge et al.ll 19771 : iFed erman et al. 1979). These theoretical (lElmegreen 



19891. Il993l : Papadopoulos et al.ll2002r) . and observational (IHonma et al.ll 19951) studies made clear that ISM 
pressure, ambient FUV field, as well as metallicity play major roles in the HI — » H2 ph ase transition. The role 



of pre ssure in particular has been highlighted over a wide range of gala xy properties (Blitz and Ro solowsky 
2004 ). and quantified in an empirical H2-pressure relation derived by Blitz and Rosolowsky ( 2006h. S uch 
a relation (herafter B-R relation) along with the well-established K-S relation kennicuttlll989[ 1998) are 



the most encompassing observational benchmarks that galaxy models must pass before they can be trusted 
in their predictions. Moreover, by linking gas and star formation (K-S relation), and the H2 phase (the 
true star formation fuel) to "macroscopic" ISM environmental parameters such as pressure (B-R relation) 
these empirical relations are natural choices for any sub-grid formulation that relates gas and star forma- 
tion in simulations of cosmological volumes where sub-grid recipes for star physics at kpc scales become 
necessary. 



Currently there is no evidence that the K-S and B-R relations hold in the extreme and very gas-rich star 
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forming galaxies discovered at high redshifts (e.g. IWalter et alj|2003l) . and there is even tentative evidence 
that the K-S relati on obtained in the local Universe may not be applicable in UV/optically selected galaxies 
at high redshifts jTacconi et alJbooj) . Detailed galaxy-sized simulations of gas and stars are thus important 
tools for exploring the robustness and possible limitations of these emperical relations in a systematic fash- 
ion. Key features of our galaxy-sized TREE/SPH numerical models of gas+stars that make them appropriate 
for such purposes are: 



• Non-equilibrium treatment of the gas thermodynamics, resulting in gas with (n, T^) ~ (0. 1 cm 3 , 10 4 K) 
(the WNM HI phase) to (n,T kin ) ~ (100cm -3 , 40 K), (i.e. the CNM HI and the resulting H 2 phase). 

• Tracking temporally and spatially varying radiation fields (profoundly influencing the HI — > H 2 phase 
transition within a galaxy) using time-dependent stellar evolution libraries. 

• A time-dependant sub-grid physical model of the HI/H2 mass exchange that readily incorporates the 
varying ambient conditions expected for the ISM within an evolving galaxy. 

• Star formation controlled by a Jeans-mass instability criterion. 



The latter is enabled by the fact that our code tracks the gravitational and thermodynamical state of 
the gas and can identify gravitationally unstable regions down to the temperatures and densities typical of 
Giant Molecular Cl ouds (GMCs). It is in such regions that strong observational evidence and theoretical 
considerations (e. g. lElmegreenll2000l. 120021) suggests that star formation occurs. Finally, note that given the 
continuous mass exchange between the WNM and the CNM gas phase, and the n qn-equihbrium cond itions 
often found for the former even for quiescent environments in the Milky Way (fWorfire et al.l 120031) . any 
successfull time-dependent treatment of the HI <-» H 2 mass exchange within evolving galaxies must track 
the ISM thermodynamics. 



2.1. The H 2 model 



The subgrid cloud structure model used by IPelupessv et all (120061) to describe the HI <-» H? m ass ex- 
change is constructed using widely observed ISM scaling laws ( Larson! 1981 : Hever and BrunJ2004 ). shown 
to hold generally for gas clouds virialized under a background pressure JElmegreenlll989h . A major devel- 
opment since our first use of this sub-grid cloud representation in o ur numerical model s is that these scaling 
laws have now been found to hold for extragalactic GMC s as well (Bolatto et al.ll2008l) . Below we describe 



its main features, while more details can be found in lPelupessy et al 



For a cloud with radius R, mean density (n), and internal density profile n(r) oc 1/r, consisting of a 
molecular core and an outer HI gas layer of a HI — > H2 transition column density A^ r (HI), under irradiation 
by an external stellar UV field the molecular fraction can be expressed as 



fn 2 



M(H 2 ) 

~m7 



= exp 



Mr ' 

'(n)R_ 



(1) 
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Here we will assume th at the gaseous ISM is composed of st ructure conforming to the widely observed 
density-size scaling relation (ILarsonlll98ll : IPelupessy et al.ll2OO60 . 



(n)R = 4.7 x 10 



21 



Pe/h 



10 4 K 



cm 



1/2 



cm 



(2) 



Hence, a calculation of the thickness of the neutral layer A^(HI) gives the local molecular fraction from 
equations Q] and |2] For this transition column density Af (HI) a differential equation can be formulated that 
describes the time evolution 



daNtr 
~dt 



= r d i s e 



-crN tI _ 



,aN a /o-N c 



1 



(3) 



where T/ = 1 / (2nR f) is the H2 formation timescale. The H2 formation rate function R f depends on temper- 
ature Tic, metallicity Z, and a normalization parameter p , (encoding the uncertainties inherent in its absolute 
value, with = no = 3.5 corresponding to the Jural (119741) formation rate of Rj ■ ~ 3 x 10~17 s" 1 at T = 100K), 

as 



R f = 3.5 x IQ-^fiZ ' h 



100 ^ 



1/2 



Sii{Tk)lH 2 cm 3 s 1 . 



(4) 



The function SniT^) expresses the HI sticking probability on a dust grain and forming an H2 molecule that 
then detaches itself from the grai n with a probabili t y 7h v Here we adopt SniT^) = [1 +(kBTf c /E a )]~ 2 (E D /kB ~ 
100 K) obtained by the study of iBuch and Zhangl (1199 lb and 7h 2 ~ 1- The dimensionless parameter 



Gk Q 
n e Rf 



(5) 



measures the relative balance of the H 2 dissociation versus the H 2 formation, with k a =4x 10 s being the 
(unshielded) H2 dissociation ra te and G the far-UV radiation field in Draine field units (2 x 10 7 photons s -1 
cm 2 between 11.2 and 13.6 e V. iDraindl 1978H . The dimensionless factor $ is an i ntegral of the self-shielding 
funct i on over the H 2 column , which encompasses the details of H 2 self-shielding (iGoldshmidt and Sternberg 



3] wi thin our dynamica l 



Pelupessv et all (12006). 



19951 : IPelupessy et a l. 2006). For a detailed explanation of the solution of Equation 
model and key dependencies of the HI — > H 2 phase transition the reader is refered to 
The numerical simulations presented in the section [3] use the solution to the fully time-dependent Equation [3] 
For the moment we will consider the equilibirum solutions first to gain some qualitative insight in the B-R 
relation. 



2.2. Equilibrium results: The B-R relation 

The equilibrium transition column density for th e fiducial case where the density, radiation field etc of 
a given patch of ISM is constant in time, is given by (IPelupessy et al.l l2006). 



v , ( 3 Gkn T 

N tr {HI) = In 1 + — — °$ 

ctfuv V 2v Rfn 



(6) 
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where v - nRapuv I i^ + nR&FUv)- Together with Equations [TJ and [2] one can calculate the corresponding 
equilibrium molecular to neutral ratio R m = fu 2 /fm = /h 2 /( 1_ /h 2 )- The latter depends only on the following 
local conditions of an ISM gas parcel: the density (n), temperature (7^), metallicity (Z) and impinging UV 
radiation field(G) as well as the local turbulent velocity field with velocity dispersion a through the total 
external pressure 

P e /k B = n (T + 54a 2 ) Kcm" 3 (7) 

(needed in Equation©. In Figure [TJwe show the resulting R m -P e relation and compare it to the observed one 
(Blitz & Rosolowsky 2006). The shaded regions indicate the scatter in the observational B-R relation, both 
for the individual measurements in a given galaxy as well as between different galaxies. For interpreting the 
plots in Figure[TJwe need to consider the following: a) the observed B-R relation is one between the projected 
molecular-atomic ratio and midplane pressure (estimated from projected quantities), and thus not exactly 
the same as the theoretical points in Figure [TJ that use the volume-averaged local R m ratio and pressure, and 
b) these results correspond to ISM equilibrium. In practice the ISM may not be even in an approximate 
equilibrium, especially for low density/low pressure regions where the timescales to equilibrium are the 
longest. For the moment we defer discussion of projection and non-equilibrium effects to the investigation 
of realistic galaxy models in Section [3] 

A number of important points becomes apparent from the panels in Figure [TJ namely: 

• for a wide range of parameters a B-R type of relation does emerge. Variations in temperature (fig. (It), 
velocity dispersion (fig. Q})), and radiation field (fig. [TJl) as well as the formation rate parameter /j, 
(fig- Et) have only a minor effect. This shows that a (B-R)-type relation is plausible from a theoret- 
ical point of view, while its various functional dependances remain within its observational scatter 
expected within and between galaxies. 

• Metallicity has a more pronounced effect on the H2-pressure relation. Figure shows that for low 
metallicity environments the theoretical relation tends to fall significantly below and outside the nom- 
inal range, while it steepens at low pressures. Systematic studies of the B-R relation in low metallicity 
galaxies may reveal such deviations. Note however that it is still possible to shift the theoretical points 
back to the nominal B-R relation by assuming e.g. a lower radiation field. 

• The plotted relations in general have slopes and normalization similar to the observed B-R relation, 
but not necessarily equal (though remaining mostly within the expected observational scatter). In- 
troducing a secondary dependence of one of the other variables on pressure can easily produce exact 
matches of the B-R slope. In a stationary model there is no unique way of doing this, though. For 
example, postulating a G dependency on pressure (G oc P 1 / 2 ) or a velocity dispersion dependency 
a oc P 1 / 4 or some suitable combination of those will result in a relation with a slope close to the 
observed value (~ 0.92). While such secondary relations are plausible (e.g. ISM environments with 
large pressures tend to host more vigorous star-formation and will thus have higher G's), it is uncertain 
whether they indeed emerge in a more realistic time-dependent galaxy-size model of gas and stars. 

• The pressure dependence of R m expressed in the observational B-R relation could effectively boil 
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down solely to a density dependence given that midplane pressures observationally are estimated 
using a constant velocity dispersion (the dominant pressure contributor in th e CNM ISM). Pr evious 



work though has deduced a direct dependence of R m on the ISM pressure (lElmegreenl Il993l) . To 



distinguish between these possibilities in Figure QJ we investigate the H2-pressure relation at constant 
density (i.e. varying the pressure only through the velocity dispersion). It can be seen that the B-R 
relation is still present but tends to flatten (especially for high densities) to oc P 5 . Additional sources 
of pressure (magnetic fields, ram pressure, shocks) may be expected to behave similarly: increasing 
R m but not as strongly as an increase in density would do it. 

Finally, unlike the investigation of B-R relation a similar one for the K-S relation must involve the full 
dynamical treatment allowed by our models given that in our approach star formation is controlled by the 
Jeans mass criterion (a dynamical one) rather than any parametric formulation. 



2.3. The CO model 



Direct detection of H2 gas is difficult given that its lowest transition S(0) : J u -Ji = 2-0 at 28/xm 
still has E2o/kB ~ 510K, much too high to be substantially excited by the typically much colder H2 gas 
(Tt ~ 15-60 K). This fact, along with its small Einstein coefficient (A20 = 2.94 x 10~ n s~ l ), diminishes its 
luminosity, while at 28/xm the Earth's atmosphere is mostly opaque, further compounding the observational 
diffuculties of its detection. These difficulties made the next most abundant molecule after H2 itself, CO 
(with [CO/H2] ~ 10~ 4 for Solar metallicities) and its easily excited rotational lines (mostly CO J=l-0 at 
1 15 GHz with Eio/ke ~ 5.5K and n K jt ~ 400 cm~ 3 ) the molecula r gas tracer of choice via the so-called CO- 



H2 conversion factor (IDickman et al.ll 19861 : iSolomon et al.ll 19871 e.g.) It must be noted that all fundamental 
relations that involve the H2 gas distribution in galaxies have been deduced for CO-bright H2 gas. This 
may not encompass the bulk of the molecular gas phase, espe c ially in metal-poor and/or FUV-aggressive 
ISM environments (|Maloney and Black! Il988l : IPak et al.lll998l : Bolatto et al.lll999i) . Such conditions can 
be found in s piral disks at large galactocentric distances because of well-known metalicity g radients (e.g . 



Henr 



V 



1998 



Madden et al. 



Gar netM Il998l and references therein), as well as in dwarf irregular galaxies (|Israellll997 
19971) . Finally metal-poor systems with significant star forma tion rates (and thus strong FUV 



radiation fields) such as Ly-break galaxies are also known at high redshifts (|Steidel et al.ll 1999n . 



Thus in our models it would be instructive to examine the specific distribution of the CO-rich, conven- 
tionally observable H2 gas. Moreover, by tracking the evolving ISM environment in which real gas clouds 
are immersed, we can identify conditions regions and epochs in which CO-bright H2 gas in galaxies may 
significantly underestimate its true distribution during the evolutionary track of a galaxy. Including it in sim- 
ulations allows a dynamic examination of the CO-H2 relation, contrasting and compleme nting those based 



on static Photo-Dissociation Regions (PDRs) models (|Pak et al.lll998 



Bolatto et al. 



1999). 



For the CO model we adopt a similar model as for H2: the size of the CO dominated region within a 
spherical clou d can be estimated by con sidering the width of the C + layer that surrounds a mixed C°, CO 
inner region (|Papadopoulos et al.ll2004l. and references therein). The dominant reaction channels for the 
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log P e>t /k b (K cm- ! ) log P e>t /k b (K cm') 

Fig. 1. — Equilibrium expectations for /h 2 //hi (points) for a range of physical parameters (section 2) versus 
the observed B-R relation and its scatter (dashed line and shaded area). In panels (a)-(e) the plus symbols 
are plotted with the same default set of parameters: T = 100K, a = 8 km/s, \i = /xq, G = 10 x Go and Z = Z . 
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formation and destruction of C + that determine its layer are: 

C° + v -> C + + <?~ 

C + + e~ -> C° + ^ (8) 
C + +// 2 -> Cff^+z/ 



Following iRollig et all (|2006r) . the radius r+ beyond which the gas is C + -dominated in spherical FUV- 
illuminated clouds with a uniform density n can be estimated from 

3 x Kr'V 1 G E 2 [i FUV A v {r + )} = n R {a c X c + 0.5 k c ) , (9) 

where X c = [C/H] = 1.4 x 10~ 4 Z, the factor £,fuv ~ 2-4 accounts for the absorption at FUV wavelengths, 
and Ej is the second order exponential integral 

roo e -y£,FuvA v {r+) 

E 2 [fruvA v (r + )]= / j dfx. (10) 

J\ I 1 

The recombination and radiative association rate coefficients for the reactions: C + +e~ — > C+v and C + +H2 — > 
CHj + v that destroy C + are a c = 3 X 10 -1 1 cm 3 s -1 and k c = 8 x 10~ 16 cm 3 s" 1 , while the extinction A v from 
r+ to the edge of the cloud R, given by (again for a n oc r -1 density profile) 

A v (r+) = 0.724 cj v Z?i/? In . (11) 

The product nR can be eliminated using the linewidth-size relation Equation |2] Equation [9] is solved nu- 
merically for r+ by simple root finding. From this we obtain the C043right part of the H2 cloud through 
fco = (r+/R) 2 . Given the fact that CO formation happens at much higher densities (n > 1000), and corre- 
spondingly shorter timescales, than H2 formation, the stationary treatment adopted for the CO chemistry is 
appropiate. 



3. The dynamical model: gas+stars galaxy simulations 
3.1. Simulation code 



The code calcu lates gravity using a TREE code (IBarnes and Hutlll986h and gas dynamics using the SPH 
formalism (see e.g. Monaghan 19921) . and the conservative formulation of ISpringel and Hernquistl ([2002). 
An advanced model for the ISM medium is used, a star formation recipe based on a Je ans mass c r iterion , 
and a well-defined fee dback prescription. The code is described and tested in detail in IPelupessyl (120051) : 
Pelupessy et all (120041) . below we will only give a brief description of the relevant physical ingredients. 



3.1.1. ISM model 



Our ISM model is similar, albeit simplified, to that of IWolfire et al.l (119951. 120031) . We solve for the 
thermal evolution of the gas including a range of collisional cooling processes, cosmic ray heating and 
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ionization. From the viewpoint of our application here the most important feature is the tracking of the 
WNM and CNM HI phases. The latter is where high densities and low temperatures allow the H2 molecules 
to form and survive, with the H2 gas phase (Section |2]) then naturally completing the ISM treatment. 

The FUV luminosities of the stellar particles, which are n eeded to calculate the photoel ectric heating 
from the local FUV field, are derived from Bruzual & Chariot (|Bruzual A. and Charlodll993l . and updated) 
population synthesis models for a Salpeter IMF with cutoffs at 0.1 M Q and 100 M Q . In the present work 
we do not account for dust extinction of UV light, other than that from the nata l cloud: for a young s tellar 
cluster we decrease the amount of UV extinction from 75% to 0% in 4 Myr (see [Parravano et al.ll2003r) . 



3.1.2. Implementation of the H2 model 

The gas particles in the code are assigned, in addition to the usual density p, internal energy u, etc, 
a varying local molecular gas fraction f m . We then use Eq. [3] to track the evolution of aN b , and thus /„, 
(through eq. [T]and eq. during a simulation timestep dt. If the temperatures are too high for H2 formation 
to occur we solve for the evolution of f m using pure photo-destruction while for T|< > 3000 K, w e treat the 



collisional destruction process of the remnant molecular gas using rates from lMartin et al.l (11998). 



The density that enters in those equations is assumed to be the mean density (n) given by the SPH 
density at the particle position, and the temperature the particle temperature (both taken constant during the 
timestep). The radiation field is calculated from the distribution of stars, assuming no other extinction apart 
from that in the natal clouds. For the macroscopic pressure P e we need the local velocity dispersion a. For 
this we take the formal SPH estimate 

°y = —to- {v)jfw(\nj\,hj) (12) 
i Pi 

with v ; and m t the particle velocities and masses, (v) ; the local bulk velocity. 



3.2. Star formation: SN feedback, and H 2 as an additional SF regulator 

The coldest and densest gas in our simulations is found in the CNM phase with the H2 formation occur- 
ing on the formation timescales of the Giant molecular clouds (GMC) complexes embedded in this phase. 
Following gravitational instabilities further "down" in the CNM phase would require additional physics (e.g. 
CO and H2O cooling of dense molecular cores, the emergence of the IMF, etc) as well as demanding much 
higher numerical resolution, currently unattainable. At this point it is appropriate to introduce a prescription 
to further track the star formation process. 

The first assumption we make is that star formation is governed by the gravitational unstability of 
gas clouds, with a region considered unstable to star formation if its local Jeans mass Mj < M rR f, wher e 



M re f S3 10 4 5 M Q is a reference cloud mass. The exact value of M re f is not important (lGerritsenlll997h . 



and provided that it is always well-resolved by our simulation this star-formation criterion precludes the 
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emergence of numerical artifacts that can result from insufficient resolution of the Jeans mass. Moreover an 
M re f smaller than typical GMC masses (as the chosen values are) makes the Jeans instability criterion select 
conditions "deep" into the dense parts of the CNM phase. This makes this criterion a good assumption for the 
onset of star formation, mirroring the irreversibility of the (gas)— >(stars) transition observed in nature once 
dense CNM clouds form. Once a region is deemed unstable it proceeds towards star formation by converting 
some fraction e$F of the gas particle to stars after a delay time. This delay is taken to be proportional to 
the local free fall time: r s f = f s ftff = The delay factor f s f is uncertain, but fro m observations a value 

? s f ~ 10 seems necessary to account for the observed inefficiency of star formation ( Zuckerman and Evans 



19741) . The actual rate of star formation is then determined by balance between gas cooling and the far-UV 



and SN heating. We will refer to this star formation model as the Simple Delay (SD) model. 

Our ISM model allows to set the local H 2 gas mass fraction as star formation regulator in the dynamical 
setting of an evolving galaxy. Irrespective whether H2 formation ahead of star formation is incidental (e.g. 
cold and dense gas forms H2 on its "way" to gravitational collapse and eventual star formation) or instru- 
mental (e.g. H2 must form first so that CO and other powerful molecular coolants can form and "drive" the 
gravitational collapse further towards denser and colder ISM regimes), this is a very important step towards 
a better, much more realistic rendering of the star formation process in numerical models. We implement 
this molecular regulated (MR) star formation by converting the molecular (f m ) mass fraction of an unstable 
(i.e. Mj < M re f) gas particle to stars (with a minimum value of f m = 0.125, corresponding to a particle mass 
of ~ 60M Q , to avoid the creation of very small star particles). Unlike the (SD) recipe that needs an adhoc 
e$F value, the (MR) one contains a physical basis for this part of the star formation modeling and thus no 
longer needs a star formation efficiency parameter e$F ■ 

The mechanical energy output of stars is reasonably well known but it has been proven difficult to 
include its feedback (i.e. supernovae and stellar winds) self-consistently in galaxy-sized ISM simulations. 
The reason for this is that the effective energy of such feedback depends sensitively on the energy radiated 
away in thin shells around the bubbles created, which would need prohibitively high resolution to follow. 
In SPH codes there have been conventionally two ways to account for feedback: by changing the thermal 
energy input and by acti ng on particle velocities. Both are unsatisfactory, as the thermal m ethod suffers from 
overcooling (|Katall992l ) and the kinetic method is too efficient in disturbing the ISM ( Navarro and White 



19931) . Here we use a method based on the creati on of pressure particles that act as norm al SPH particles 



in the limit of the particle mass m p going to zero (|Pelupessy et al.ll2004l : IPelupessyl 120051) . Such a pressure 
particles is associated with every newly formed star particle and will receive its feedback energy, acting on 
the surrounding gas particles through the usual SPH particle forces in the limit that m p — > while simul- 
taneously keeping the product of particle mass and specific thermal energy, m p x u p , fixed. The thermal 
evolution (the time dependence of m p x u p ) is specified by adiabatic expansion and the energy input from 
young stars. For this energy injection rate we take E = e sn n sn E sn /At, with E sn = 10 51 er g, e Jvj — •!> n sn — 

0.009 

per M Q and At = 3 x 10 7 yr. The efficiency e sn thus assumes that 90% of the energy is radiated away in 
thin, dense shells. 
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3.3. Galaxy models 



The galax y models we use stem from the analytic disc galaxy models of iMo et all (119980 . constructed 



as described in Spring el et all (|2005). They consist of a disk consisting of a stellar and a gaseous component 



embedded in a dark halo. The stellar disk has an exponential disk radial profile (with scale length R+), 

PdM(R,z) = ^ exp(-R/RJsech 2 (z/h z ). (13) 

The gas disk is set up in vertical hydrostatic equilibrium with a surface density profile consisting of an 
exponential component (so proportional to the stellar density) and a more extended component 

S = S g /(l+R/R g ) (14) 

cutoff at a radius 8 x R g . Note that this distribution is necessary to match closer the observed gas distributions 
that typically extent well beyond the stellar ones. Apart from the radial profiles the stellar and gas disk are 
initialized as smooth initial conditions and the gas is setup with a constant temperature (8000K). Finally the 
dark halo has an Hernquist profile 

Mhalo a 

Phalo(r)= — 3 . (15) 

The Hernquist scale parameter is related to the more f amiliar scale param eter r s and the concentration index 



c of Navarro, Frenk and White (NFW) profiles (see ISpringel et al.H2005l . for details). We do not include a 
bulge component here. 

We take models of different size by choosing the total baryonic mass, the disk being a mass fraction 
/baiyon of the total mass, and consider galaxies ranging in mass from Mbaryon = 10 8 M Q to Mb ary0 n = 10 9 M Q , 
with /baryon = 0.041. The smaller scale is representative of dwarf galaxies and the bigger of small disk 
galaxies. We vary the gas mass fraction of the disk from / gas = 0.5 (for the low mass model) to / gas = 0.2 (for 
the high mass model), roughly mirroring the observed correlation between gas fraction and galaxy size. The 
total halo virial mass is fixed by M vu = Mbaryon //rf (M v i r = M na io +M\, aryon ), which gives the virial velocity 
and radius through the relations 

v 3 . 

M vir = v -^, (16) 

l0GH(z) 

R ™ = ioGH( Z y (17) 

(18) 

assuming a virial overdensity A = 200. Determining the c index gives the halo scale length. The metallicity 
of each model will be taken to be constant during the run, but we will consider different metallicities, namely 
models at solar metallicity Z Q and at 0.2 x Z Q . These models are summarized in Table[T]as models Al to 
CI. Other parameters of the models are not varied here: we fix the spin parameter at A = 0.05. The A 
implicitly determines the scale lengths of the gas and stellar disk, as the angular momentum in the disk is 
assumed to scale with the angular momentum in the parent halo. The scale height of the stellar component 
is taken to be a fraction of the radial scale length: Zh = 0.2-0.4 x R+. The galaxy models are realized with 
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mass resolutions for the gas particles of 200- 10 3 M , also given in Table [T] Finally two models with very 
high gas mass fractions are also run (Dl and El), representing extreme systems expected at early epochs of 
galaxy evolution. For these we set an initial gas mass fraction of / gas = 0.99. 

3.4. Runs 

In addition to different galaxy models A 1 -El we also test our two different star formation recipes (MR 
and SD) described in Section I3T21 Each model is run well beyond the time strong evolutionary effects take 
place (investigated in section @]) and until a dynamic equilibrium for the star formation sets in. At this 
evolutionary point, i.e. ~ 1 Gyr of simulation time after the start of the simulation, we analyze the resulting 
gas distributions. Given that at present chemical enrichment effects (influencing ISM thermodynamics and 
H2 formation) and cosmological infall are not included, evolving our models for much longer timescales is 
of limited value, simply resulting in a steady depletion of their gas reservoirs. 

3.5. Results 

Table|2]gives an overview of the molecular gas fractions and star formation rates of our runs. From this 
table it can be seen that the low metallicity models (Z = Z Q /5, A B and E models) have a molecular fraction 
f m ~ 0.03-0.05, while the models at solar metallicity (C and D models) reach up to f m ~ 0.4-0.6. Models 
with the same metallicity reach similar molecular fractions while the structural parameters of the galaxy 
models seem to have only a minor influence on the global molecular fraction, at least over the limited range 
explored here. It is also important to point out that the SPH particle mass has little effect on the basic physical 
quantities examined here, suggesting that adequate numerical resolution has been reached to describe the 
physical mechanisms considered. 

A comparison of SD and MR star formation recipes shows that the latter results in less molecular gas 
and increased SFRs. A useful measure of the star formation efficiency is the gas consumption timescale r x = 
M X /SFR, also given in Table 12 and calculated separately for atomic and molecular hydrogen. For the low 
metallicity models we find typically thi ~ 5-8 Gyr and t# 2 ~ 0.15-0.3 Gyr, while for the high metallicity 
models the timescales for HI and H2 become comparable, with tri ~ 1.5-3 Gyr and t# 2 ~ 0.6-5 Gyr. These 
results seem largely independent of the SF model adopted and reflect a general characteristic of low-Z versus 
high-Z systems, namely that the former are much more WNM-dominated than the latter. This can be seen in 
the last column of Tabled where the Mcnm/Mwnm ratio is also tabulated (CNM is taken to be all gas colder 
than 1000 K, WNM gas with 1000 < T < 10000). For systems that are WNM-dominated, hydrogen will be 
overwhelmingly atomic, and the large disparity between the HI and H2 gas consumption timescales simply 
reflects the one between atomic and molecular gas reservoirs and the fact no Jeans-unstable regions occur 
in the WNM phase and thus star formation can never directly "consume" its mass. For high metalicities the 
mass allocation between WNM HI and CNM HI and H2 becomes more even and so are the corresponding 
consumption timescales. Note that for the SD model alone one might be tempted to conclude that the short 
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H2 gas consumption timescale at low Z (or equivalently the low f m ) has something to do with the fact 
this model is formulated independent of molecular gas (thus unrealistically converting gas into stars before 
molecules can form). This is not the case: the MR model has lower consumption timescales at low Z (partly 
due to a general higher SF), and lower molecular fractions. It seems that even in the MR model the outcome 
of the star formation model is not constrained by the chemistry of H 2 formation but by the conversion of 
WNM to CNM. The difference between high metallicity gas and low metallicity gas is that molecular gas 
forms in smaller reservoir of CNM gas. Once there, evolution to star formation occurs faster than the large 
scale processes driving gas down to the CNM, so it is not necessary for a large reservoir of H2 to form. 

Compared with the SD model the MR star formation model has a ~ 50% lower molecular fraction 
and a ~ 50% higher rate of SF. This increase in efficiency of SF in terms of its molecular mass means that 
the H2 gas consumption timescale is a factor ~ 3 - 4 times shorter. The reason for the smaller amount of 
H2 gas is that the MR recipe selects SF sites that are on average denser CNM regions (where H2 forms), 
and that these regions are then converted to stars, resulting in a more efficient consumption of molecular 
gas. At least globally, M(HI+H2)/SFR seems relatively insensitive to the SF recipe chosen. Given that stars 
form unequivocally only out of the H2-rich regions of the CNM phase this suggests that any self-regulating 
mechanism responsible for distributing the gas between the SF and the non-SF phase remains broadly similar 
in these two SF recipes. Of course we must reiterate that the MR star formation recipe is the more physical 
of the two, doing away with the adhoc £sf parameter typically used in numerical models. The emergence 
of a robust global efficiency M(HI+H2)/SFR out of dynamic galaxy models where only cold, dense, and 
H2-rich gas is allowed to form stars confirms the trustworthiness of (K-S)-type phenomenological relations. 
This is because the latter often relate the total gas mass to star formation, irrespective of its thermodynamic 
state or molecular gas fraction. 

In Figures [2] and [3] we show the gas distribution from the simulation snapshots of models Al-El. The 
top row of Figure[2]shows the HI gas maps obtained from projecting the neutral mass fraction of the particles 
for runs using the SD star formation model. In the middel row the molecular gas distribution obtained from 
f m directly is shown, while the bottom row shows the maps for the CO-rich H2 distribution (determined 
as in Section |2~3T ). Figure [3] shows the equivalent maps for the MR star formation runs. Some features of 
these maps are common across different models. For example the panels for the CI run (spiral galaxy/solar 
metallicity, middle panels) show the frothy appearance of the neutral gas distribution typical for a star 
forming ISM. Comparing the HI and H2 maps we see from the enhanced contrast of the molecular map that 
the H2 tends to concentrate in the higher density regions. In the outer galaxy regions the H2 distribution cuts 
off before the HI distribution and the smoothness of the gas distribution shows little feedback from stars. 
The third row panels shows the CO distribution concentrated towards high column densities in the central 
regions and dense clumps. Much the same pattern is visible for the equivalent low metallicity run Bl. Low 
metallicity means less H2 and CO formed, with CO restricted to the very highest density clumps. Note also 
that there is a big smooth region where star formation and H2 are absent, a pattern repeated in the Al run 
(note that the linear scale of the maps are different). The pure gas models (two right most panel s in figs. [HO 
are stable despite the high gas content, stabilized by supernova feedback (|Springel et al.l l2005). 

The resulting gas distribution for the MR star formation model (fig. [3]) is generally similar to the SD 
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Table 1 : Overview of galaxy model parameters. The letters A-E indicate different structural properties and/or 
metallicities while numbering indicates different resolutions used on otherwise identical models. The gas 
distributions of models Dl and El consist of equal mass exponential and extended disk, the other models 
consist purely of the extended distribution eq. [14] 



Model 


^baryon 




/gas 




Z 


Al 


10 8 


2.3 x 10 9 


0.5 


200 


0.2 


A2 


10 8 


2.3 x 10 9 


0.5 


500 


0.2 


Bl 


10 9 


2.3 x 10 10 


0.2 


500 


0.2 


B2 


10 9 


2.3 x 10 10 


0.2 


1000 


0.2 


CI 


10 9 


2.3 x 10 10 


0.2 


500 


1. 


Dl 


10 9 


2.3 x 10 11 


.99 


1000 


1. 


El 


10 9 


2.3 x 10 11 


.99 


1000 


0.2 



Table 2: Results: molecular masses and fractions for runs. 



Run 


SF 


M gas 


M m 


Mh 2 


/h 2 


R 


SFR 


THI 


T~H 2 


JCNM 
fwNM 




model 


(1O 7 M ) 


(1O 7 M ) 


(1O 7 M ) 






(M /yr) 


(10 9 yr) 


(10 9 yr) 




Al 


SD 


4.5 


4.0 


0.2 


0.045 


0.05 


0.005 


5.8 


0.3 


0.64 




MR 


4.4 


4.1 


0.1 


0.023 


0.024 


0.006 


4.7 


0.11 


0.41 


A2 


SD 


4.4 


4.0 


0.21 


0.047 


0.052 


0.005 


6.0 


0.31 


0.64 




MR 


4.2 


3.8 


0.1 


0.024 


0.027 


0.006 


4.7 


0.13 


0.4 


Bl 


SD 


19.0 


17. 


0.86 


0.045 


0.049 


0.01 


11.8 


0.58 


0.61 




MR 


19.0 


17. 


0.4 


0.02 


0.022 


0.015 


7.8 


0.17 


0.39 


B2 


SD 


19.1 


17. 


0.84 


0.044 


0.048 


0.009 


12.3 


0.59 


0.59 




MR 


18.6 


17. 


0.34 


0.02 


0.02 


0.015 


8.4 


0.16 


0.38 


CI 


SD 


18.6 


6.7 


11. 


0.61 


1.7 


0.014 


3.2 


5.4 


5.28 




MR 


17.1 


11. 


5.4 


0.3 


0.49 


0.031 


2.5 


1.2 


1.81 


Dl 


SD 


88. 


41. 


46. 


0.52 


1.1 


0.14 


2.0 


2.2 


3.06 




MR 


76. 


53. 


20. 


0.27 


0.4 


0.24 


1.6 


0.6 


1.32 


El 


SD 


90. 


83. 


4. 


0.045 


0.05 


0.1 


4.9 


0.24 


0.59 




MR 


89. 


84. 


1.8 


0.02 


0.02 


0.17 


3.5 


0.08 


0.34 
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model, especially on large scales. On small scales the H2-regulated SF structures are affected by feedback 
due to a star formation more biased towards high density regions, which results in a more bursty star forma- 
tion mode with large "bubbles" (like that seen in the centre of the panel of El model) forming more often. 



3.5.1. N on- equilibrium HI <-> H2 mass exchange 

For typical CNM HI densities the H2 formation timescales are ~ 5 - 50 Myr (Eq. 4). The frothy gas 
disk structures evident in Figures [2] and [3] indicates the dynamic nature of the processes. In order to show 
this it is necessary to resolve the feedback effects in enough detail. Together with an assumed equilibrium 
f m this will result in a more static picture for the ISM. 

In figure@]we demonstrate the non-equilibrium between the HI and H2 gas phases by plotting the equi- 
librium f m from section |Z21 as a function of density, and the actual non-equilibrium f m from the simulations 
(for Z Q and Z Q /5). As we can see from this snapshot the gas is out of HI/H2 equilibrium for a large frac- 
tion of gas particles, with only the highest densities converging to equilibrium. The "ergodic" rather than 
particle-ensemble representation of this is demonstrated by the evolutionary track of a single gas particle 
in the f m -n plot (dashed line) showing that it spends most of the time away from equilibrium areas. This 
also means that both the collapse to higher densities and the destruction of H2 in the diffuse gas phase is 
slow compared to the evolution of gas structures and ambient ISM conditions (further justifying the use of 
a time-dependent computation of f m ). 



3.5.2. The CO tracer molecule versus H2 



The CO-rich H2 distribution in Figures [2] and [3] is markedly different from the total H2 distribution 
demostrating a variable CO-to-H2 mass ratio, especially for metal-poor systems. This can be seen more 
clearly in Figure [5j where we plot the pixel values drawn from the H2 distribution against those drawn 
from the CO map. The values in figure [5] are scaled so that the maximum pixel values of the CO and the 
H2 distributions match, which could be considered as an effective calibration of the CO-to-H2 conversion 
factor for our simulations (observationally also performed at the bright end of CO-luminous H2 clouds). It 
is clear that the CO-to-H2 factor is not constant: the CO-rich H2 distribution does show a tight correlation 
to the total H2 mass distribution, but one that is much steeper than linear (and thus difficult to calibrate 
observationall y without running into sensitivity limitatio ns). This has been suspected and argued widely in 
the literature (Ma loney and Black!ll988l : |Pak et al.N 19981. e.g.), so it is interesting that in our simple model 
this is demonstrated in the context of evolving galaxy models. Note that although a steeper CO-H2 relation 
seems to appear for metal-rich compared to the metal-poor systems for the SD simulations (contrary to 
what would be expected for stationary cloud/radiation field models, e.g. Maloney and Black1ll988 : Pak et al. 
1998), this trend dissapears in the more realistic MR models. 



The non-linear relation between CO and H2 may also raise the possibility that (especially metal poor) 
galaxies can be "CO-dark" during certain epochs of their evolution (e.g. immediately after a burst of star 
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Fig. 3. — Gas distribution of simulations. Shown are (from top to bottom) HI distribution, the H2 distribution 
and the CO map for runs (from left to right) Al-El after 1 Gyr of evolution and MR star formation. 
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formation and the subsequent enhancement of far-UV radiation), while H2 gas is still there continuing 
forming stars. Such systems may appear as having much larger than usual star formation efficiencies, i.e. 
forming stars out of seemingly very little CO-bright H2 gas. (we return to to this point in section U.5 .3 1 and 
S. 

Even strongly varying CO-to-H2 relations can be difficult to discern and calibrate observationally. Di- 
rect methods would entail independent observations of CO and H2 at comparable resolutions, an improbable 
proposition given that direct H2 observations are difficult and thus rare. The latter are: a) observations of 
its l owest excitation S(0) line (IValentijn and van der Werllll999l. e.g.) which can be excited for CO-deficient 
H2 (IPapadopoulos et al.ll2002h . b) H2 absorption studies in the far-UV. The S(0) line emission observations 
at 28^m are restricted to Space, where the small apertures deployed until now cannot match the resolution 
or gas-mass sensitivity attainable with ground-based CO observations using mm/sub-mm telescopes. Line 
absorption studies on the other hand are restricted by nature to single and special lines of sig ht, making rou 



tine c omparisons with CO observations difficult especially for extragalactic environments. iKrumholz et al. 



(|2009r) find considerably more H2 detected by far-UV absorption studies with FUSE at low column densi- 
ties than their CO, H2 equilibrium model predicts, indicating that a CO-deficient diffuse H2 phase is indeed 
present in the Galaxy. An indirect, but nevertheless powerful method relies on CO and C + observations at 
158/im where any "excess" C + emission, after correcting for contributions from H + and WNM, CNM HI gas 
phases, is attributed to H2. Such observations have indicated ~ 10- 100 times more H2 gas than what CO 
emission reveals in the metal-poor and far-UV intense environment of the dwarf irregular IC 10 (Madden 
et al. 1997). These are indeed the type of environ ments where the lar gest disparities between H2 and CO- 
rich H2 distributions are expected from static (e.g. iBolatto et al.lll999h as well as our own dynamic models 
(see Figures 2, 3). Neverthless C + observations still suffer from similar limitations like those of the S(0) 
line with the aforementioned example being one of the very few cases of meaningful comparisons with CO 
observations. Finally, with the CO-rich regions restricted deeper and deeper into H2 clouds as metalicities 
decrease (Figures 2, 3) a rising "overpressure" on these regions is ex pected by the overlying CO-deficient 
H2 gas. Such an effect has been recently detected dBolatto et al. 2008 ). 



3.5.3. The K-S empirical relation 



Past investigations of the K-S em pirical relations in galaxies used only stationary models (|Dopita and Ryder 



1994 : iRobertson and Kravtsovll2008l) . Testing for the emergence of robust K-S relations using dynamic 
galaxy models like ours could thus be very interesting, especially when the much more realistic H2-regulated 
star formation recipe is implemented. The relations between the gas mass surface densities (total, HI, H2 
and CO-rich H2) and the local star formation density in our simulations are shown in Figures [6] (the SD 
model) and |7] (the MR model) (th e CO s urface density is normalized as in Figure [5). In each case the K-S 
relation as expressed in IKennicuttl (11998!) is also shown. 



The first notable characteristic is that the K-S relation for the total gas or the HI surface density is 
closer to the observationally derived one for all the systems simulated here, both in terms of normalization 
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and slope. The slopes are generally between n m 1 .4 and n « 2. In addition, the K-S relations for the CO- 
rich H2 gas (tracing the densest molecular phase present in our simulations) have a more shallow slope (n rj 
1-1 .4). This tendency of the K-S relation towards more linear slopes as one progresse s from HI to CO-bright 
and then HCN-bright H 2 (n > 10 5 cm~ 3 ) gas has also been noted observationally dwong and Blitz 20021 : 



Gao and Solomon 



2004). The CO molecule will form in the highest density peaks of the H 2 distribution 



where short dynamical scales make them most intimately linked to the star formation. 

There is some remaining uncertainty in the normalization of our derived empirical K-S relations when 
it comes to those involving the CO-rich H2 gas (and thus directly comparable to observations). This is 
because far-UV absorption from the intervening dusty ISM (besides that in the natal clouds we considered 
here) will affect CO much more than the self-shielding H2 throught the galaxy, making the CO-rich H2 
distribution more extended than depicted in our models. This is expected to "shift" all our derived K-S 
relations involving CO to the right, but it is unlikely to reduce significantly the large deviations we find for 
the K-S relation in our simulated metal-poor systems (Z = O.2Z ) where far-UV is much less absorbed in 
their intervening ISM (i.e. between molecular clouds). In such systems the corresponding K-S relation is 
shifted upwards by a factor of ~ 5 - 10 with respect to the metal-rich ones (Fig 6, 7), yielding a much more 
efficient star formation per CO-rich H2 mass. This is because in metal-poor systems H2 and CO can form 
only in much denser gas (where star formation is most efficient) deeper in the CNM clouds, where higher 
densities make up for the loss of dust surface for H 2 and eventually CO formation. A similar effect was 
noted in metal poor dwarf galaxy IC10 (ILeroy et al.l l2006). 

It is notable that the high gas fraction models (Dl and El) have a somewhat lower SFR, i.e. their star 
formation rate per unit surface at a given surface density in these systems is lower, although the difference 
is not large (< 50%). This is likely due to the fact that between two systems with similar total surface 
densities, the gas-rich one would necessarily have a smaller stellar component constraining the gas in the 
z-direction than the gas-poor one, resulting to less star formation per total surface density. Such effects 
have been d escribed in the past, us ing K-S relations that involve the gaseous as well as the stellar mass 
component ( Dopita and Ryder 1994J). Finally, the MR models show a slight upward shift of all the K-S 
relations, but otherwise similar behaviour, despite the fact that the star formation is formulated in terms of 
the local molecular gas fraction. A K-S like law is known to arise under a wide set of conditions when the star 
forma tion rate scales with the reciprocal of the dynamical timescale 1 / y/AitGp dSchaye and Dalla Vecchia 



2008). However the MR model inserts an additional non-trivial criterion (the H2 richness of the star forming 



gas), and thus there is no reason for expecting so similar results (we return to this point in the Discussion). 



3.5.4. The ^-pressure relation 

In Figure[8]we show the Hi-pressur e relation for models Al-El, where the average R m is plotted versus 
midplane pressure P ext , as estimated by [Blitz and Rosolowskyl (120060 (their equation 5), 



P ext = 272 cm" 3 K 



M pc" 



M pc" 



0.5 
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The v elocity dispersion a and the stellar scale height h* are taken to be constants (as in Bli tz and Rosolowsky 



2006), while the gas surface density and stellar surface density are derived from the projected gas and stellar 
distributions. The R m plotted is derived either from H2 or the CO-rich H2 surface density. As discussed 
previously the CO-rich H2 surface density has an uncertainty in the absolute scaling, translating in some 
arbitrariness in the vertical scaling of the corresponding R m -P ext relation, but the trend of the R m dependence 
on P ext is not affected by this. Compared with the equilibrium H2-pressure relation of Figure Q] the plots 
here are more comparable to the actual observed quantities, given that the pressure used is the same indirect 
pressure estimate used in observations. 

In Figure [S^-d a pressure dependence much flatter than the observed ^-pressure relation is found. For 
solar metallicity the trend with pressure is very flat while a t lower metallicity it is steeper, but still short 



of a n = 0.92 slope derived by blitz and Rosolowsky (2006). One the other hand, the relation derived for 



the CO-bright H2 mass fraction Rco has a steeper dependence and is very close to the observed slope for 
Zq , while somewhat steeper for the Z / 5 simulations. From what we have seen in figure [5j this is to be 
expected: this difference between H2 and CO-rich H2 in the H2-pressure relation derives from the bias of 
CO to form at higher gas densities, and thus pressures. The steeper fall-off of the CO-rich H2 mass at lower 
densities compared to that of H2 translates to a steeper fall-off of Rco at lower pressures. 

Our results indicate that, like the K-S relation for the CO-rich H2 phase, the CO-rich ^-pressure re- 
lation should show a s t rong dependence on metallicity, and for Z = Z /5 it is considerably steeper than the 
Blitz and Rosolowsky (feoOot) relation. Although hampered by the very small numb er of actual CO detec 



tions at low metallicities, the available data suggest n o strong trends with m etallicity (Bli tz and Rosolowsky 



2006). Although some shift downwards i s expected (IKrumholz et al.ll2009h . the slope should be similar. For 



example, the IC 10 data points shown by Bli tz and Rosolowsky (2006) lie very close to the normal relation. 



It will be interesting to see whether this will be borne out by further examination of low metallicity dwarf 
galaxies. In this context it is worth pointing out that the pressure estimate, based on the equation. [191 may 
not be valid for such systems due to their low stellar surface densities. 



4. Discussion 

Our numerical models suggest that after a dynamic equilibrium sets in the differences between MR 
and SD star formation recipes shown by the simulations are small. This is not only apparent from the star 
formation rates and/or a cursory examination of the resulting gas morphologies - which are determined 
mostly through feedback processes - but also borne out by the quantitative aspects of the SF behavior as 
revealed by the emergent K-S and B-R relations in our models. For the star formation process, within the 
range of SFRs explored in our simulations, the following picture then emerges: for diffuse gas to transform 
into star-forming gas it must evolve to higher densities and the conversion into H2 proceeds passively at 
a certain high density threshold (set by formation rate and metallicity). In other words: star formation is 
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biased towards the higher densities where H2 will form fast, and this becomes even more pronounced in 
low Z cases. However in the MR star formation model the phase converted to molecular gas has a star 
formation rate sufficiently slow so that its conversion into molecular gas is not the SF bottleneck (and thus 
star formation will not be affected, even in the cases where H2 formation proceeds slowly). 

We have demonstrated that robust phenomenological relations linking gas content to star formation (K- 
S relation), and molecular gas fraction (the SF fuel) to ambient ISM pressure (B-R relation), akin to those 
found in the local Universe, emerge out of realistic settings of dynamically evolving galaxy-sized systems of 
gas and stars. However, before considering the hereby presented investigation as validating the widespread 
practice of using K-S type of relations or SD star formation models as sub-grid physics in large scale galaxy 
evolution and cosmological models, we must point o ut that a considerable fraction of the c urrent stellar mass 
has been asssembled in major ULIRG type mergers (|Smail et al 1 ll997l : iHughes et al 1 ll998h . In such systems 



destruction of H2 can have much shorter time scales, induced by faster variations of the ambient far-UV field 



expected in starbursts (Par ravano et al.l l2003) yet also competing against a faster H2 formation in their high 
density ISM gas. In such settings large deviations from (K-S)-type relations could occur, given that the latter 
seem to robustly emerge only after a dynamic equilibrium between HI, H2 gas phases and stars has been 
established. Finally, in most numerical models the sub-grid K-S relation is set to use a much warmer gas 
phase (~ 10 3 - 10 4 K) thermodynamically far removed from the CNM HI and the H2 gas that are directly 
linked to star formation. 

Modelling of ULIRG type systems is computationally more challenging as the high temperatures at 
high densities mean that timestepping of the simulation will be slower and the implicit assumptions that we 
make about the transperancy of the interstellar medium break down as the star formation sites in these are 
obscured on galactic scales, so at the moment directly testing this is not possible. However, as an instructive 
first step to examine the validity of the K-S relation during the early evolution of galaxies, we can examine 
our galaxy models during the early phases of the simulation, when the ISM phases and star formation 
have yet to establish a dynamical equilibrium. This is shown in in Figure |9j from which it can readily be 
discerned that at early evolutionary timescales significant deviations of the star formation from that expected 
from the K-S relation do occur, and are especially pronounced for metal-poor systems (Al, Bl, and El). In 
the latter cases there are periods when the CO-derived K-S relation will overestimate or underestimate the 
underlying star formation which, for gas-rich and metal-poor systems (El), can last well into later evolution 
times (T ~ 0.2-0.5Gyrs). This seems to be an effect of the greater sensitivity of CO destruction at low 
metallicities where this molecule survives only in the densest of the CNM gas, itself spawing star-forming 
regions very fast, which in turn destroy CO. The fact that this behavior emerges for both a small (Al) and 
a 10 x larger metal-poor system (Bl) suggests that this "oscillating" of the SFR with respect to a K-S(CO) 
relation (Figure 9) is not due to large stochastic scatter from a smaller number of star forming sites. On the 
other hand the CO-derived K-S relation seems to remain a good predictor of the underlying star formation, 
even at early evolutionary times, for metal-rich systems with moderate amounts of gas (CI, CI -MR), i.e. 
like those used for its establishment in the local Universe. 

The largest deviations between the SFR predicted from the K-S relation and the actual one occur for the 
K-S(HI+H2) and K-S(H2) relations during early evolutionary times (T < 0.2-0.3Gyrs), and are particularly 
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pronounced for the very gas-rich sustems (Fig.9: El, Dl models). In these cases even the K-S(CO) relation 
underpredicts the true SFR, even for the metal-rich system (Dl) where CO tracks H2 well, and thus cannot 
be due to CO failing to trace H2 (the SF fuel) well. During those early epochs gas-rich systems can appear 
as undergoing periods of very efficient star formation (i.e. little CO-bright H2 gas but lots of ongoing star 
formation), where application of the K-S(CO) relation using their observed SFRs would imply much more 
molecular gas than there is. Such systems, though mo re massive that those modeled here, may have been 
recently observed at high redshifts (|Tacconi et al.ll2008r) . 

The failure of the standard K-S(HI+H2) relation to track SFRs during the early evolution of very gas- 
rich systems is rather expected given that this relation is "blind" to the thermodynamic state of the gas, 
and thus can equally well make stars out of 10 4 K WNM HI or 30 K H2 gas. Only at later times its SFR 
predictions become valid, result of a dynamic equilibrium among the various ISM phases and the stellar 
component being established. This could have implications for modelling of the gas-rich galaxies found in 
the distant Universe, or systems where major gas mass accretion events "reset" their evolutionary states back 
to gas-rich ones. In such cases the non-equilibrium, non-linear, mass/energy exchange between the various 
ISM phases and the stellar component may come to dominate significant periods of intense star formation 
and stellar mass built-up during which not even the most realistic, CO-derived, K-S relation seems applicable 
(Figure |9j models El, Dl). 



5. Conclusions 

We use our time-varying, galaxy-sized, numerical models of gas+stars that track the ISM thermody- 
namics and the HI <-> H2 gas phase exchange, to investigate: a) the emergence of two prominent empirical 
relations deduced for galaxies in the local Universe: the Kennicutt-Schmidt (K-S) relation and the H2- 
pressure relation, b) the effects of a more realistic H2-regulated star formation recipe, and c) the evolution 
of very gas-rich systems. Our models now include a separate treatment for formation and destruction of the 
H2-tracing CO molecule, which allows a direct comparison of such models with observations, and a new 
independent investigation of the CO-H2 concomitance in the ISM of evolving galaxies. Our findings can be 
summarized as follows 

• For ISM states of H2/HI equilibrium, an ^-pressure relation close to the one observed robustly 
emerges for a wide range of parameters, with a strong dependance mostly on metallicity. For the 
more realistic non-equilibrium H2/HI states only the CO-bright H2 phase shows an ^-pressure rela- 
tion similar to the one observed. 

• The H2-regulated star formation model successfully models star formation without the adhoc param- 
eter of the local star formation efficiency adopted by most galaxy-sized numerical models, while 
incorporating a fundamental aspect of the star formation process. 

• A comparison between numerical models using the usual simple-delay (SD) and the new molecular- 
regulated (MR) star formation recipes reveals very few differences. It shows a factor of ~ 3 - 4 more 
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efficient star formation per CNM gas mass than the case of MR star formation. 

• We find little sensitivity of the global SF efficiency M(HI+H2)/SFR to the SF recipe chosen, once 
dynamic equilibrium between ISM phases and stars is established, yielding confidence to (K-S)-type 
of relations emerging as a general characteristic of galaxies. 

• A non-equilibrium HI <-> H2 gas mass exchange is revealed taking place under typical ISM conditions, 
demonstrating the need for a full dynamic rather than stationary treatment of these ISM phases. 

• The CO molecule can be a poor, non-linear, tracer of the true underlying H2 gas distribution, especially 
in metal-poor systems, and even in those with very high gas mass fractions (more typically found at 
high redshifts). 

• A K-S relation robustly emerges from our time-dependent models, irrespective of the SF recipe used, 
after a dynamical equilibrium is established (T >1 Gyr). The CO-derived K-S relation has a more 
shallow slope than the one involving the total gas mass, and as in the ^-pressure relation, a strong 
dependance on metallicity is found. 

• At early evolutionary timescales (T < 0.4 Gyr) our models show significant and systematic deviations 
of the true star formation from that expected from the K-S relation, which seem especially pronounced 
and prolonged for metal-poor systems. These deviations occur even for the CO-derived K-S relation 
(the more realistic one since CO is directly observable and traces the densest H2 gas which "fuels" 
star formation), and even for metal-rich systems where CO tracks the H2 gas well. 

• The largest deviations from the K-S relation occur at the earliest evolutionary stages of the systems 
modeled here (T < 0.2Gyr) and for the most gas-rich ones. During this time significantly higher 
star formation rates per CO-bright H2 gas mass occur, and such star-forming galaxies may have been 
already observed at high redshifts. 

Finally we must note that when it comes to the gas-rich galaxies accessible to current observational 
capabilities at high redshifts, our results, drawn for much less massive systems, remain provisional. Nev- 
ertheless for more massive gas-rich systems the larger amplitudes of ISM equilibrium-perturbing agents 
(e.g. SNs, far-UV radiation fields), and the shorter timescales that will characterize their variations are more 
likely than not to exaggerate the deviations of true star formation versus the one derived from (K-S)-type 
phenomenological relations. A dedicated observational effort to study such galaxies at high redshifts (soon 
to be dramatically enhanced by ALMA), as well as extending detailed numerical modeling of gas and stars 
to larger systems (as computational capabilities improve), can help establish whether (K-S)-type relations 
remain valid during most of the stellar mass built-up in galaxies, or only emerge after dynamic equilibrium 
has been reached during much latter evolutionary stages. 
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Fig. 4. — Equilibrium and non-equilibrium f m values versus gas density. The equilibrium f m of particles 
in the simulation (diamonds) and the non-equilibrium f m (dots) as a function of density for Z = 0/5 (left 
panel) and Z = Z Q (right panel). The equilibrium f m is calculated from the instantaneous particle properties. 
The dashed line marks a time track of a typical gas particle experiencing a cycle of collapse, star formation 
and reheating, with square symbols placed at 0.93 Myr intervals. 
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Fig. 5. — H2 versus CO-rich H2 mass surface density. Plotted are the pixel values of the projections of the 
gas distribution for the SD SF model for low (left upper panel) and solar (right upper panel) metallicities 
(models Bl and CI). Lower left and Lower right panels show the corresponding results for the MR SF 
recipe. The CO pixel values are normalized as described in the text. For reference a long dashed line with 
linear slope (n = 1) is plotted, as well as a short dashed line with a slope n = 2. 
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Fig. 6. — Kennicut-Schmidt laws: The star formation density for the simulations Al-El (from top row to 
bottom row) using the SD star formation recipe vs (panels from left to right) the total, the neutral, H2 or 
CO-rich H2 gas surface density. Shown are mean relations (points are averag ed in equal logarithmic bins of 
surface density). The short dashed line shows the empirical iKennicuttl (U998) relation with a slope of n = 1.4 
as a reference, whereas the long dashed and dotted lines have of n = 2 and n = 1 respectively. Note that some 
outliers are probably affected by image artefacts from the projection of the star particles. 
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Fig. 7. — Kennicut Schmidt laws. Same as figure [6l but for the MR star formation recipe. 
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Fig. 8. — Simulation results for R m = /h 2 //hi- Panels show fi-masa function of midplane pressure for the 
low metallicity models (A1,B1 and El) for the SD star formation model (panel a) and the MR model (b) and 
for solar metallicity models (CI and Dl) with the SD model (c) and MR model (d). The star like symbols 
indicate that R m is estimated from the H2 maps, while open symbols indicate that CO is used. 
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Fig. 9. — Time dependence of the SFR versus that expected from the K-S relation and the instantanous gas 
content. Shown are the actual total star formation (drawn black line) taking place in the models Al, Bl, CI 
for SD or MR star formation recipe (indicated in the lower right of each panel). All other lines give the SFRs 
expected from various applications of the K-S relation (with slope n=1.4) using: total gas surface density 
(dashed), H2 gas surface density (dash-dotted), and CO-bright H2 gas surface density (dotted). The KS star 
formation lines are normalized such that their average after 300 Myr matches the average star formation 
over the same period. 
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